library(tidyverse)
library(rlang)
library(leaflet)
hospitals <- read_csv("../data/hospital_list.csv")
Parsed with column specification:
cols(
.default = col_character(),
zip = [32mcol_double()[39m,
`Phone Number` = [32mcol_double()[39m,
score = [32mcol_double()[39m,
longitude = [32mcol_double()[39m,
latitude = [32mcol_double()[39m,
population = [32mcol_double()[39m,
poverty = [32mcol_double()[39m,
percent = [32mcol_double()[39m,
poverty_percent = [32mcol_double()[39m
)
See spec(...) for full column specifications.
hospitals
hospital_locations <- hospitals %>%
select(state, longitude, latitude)
write_rds(hospital_locations, "../data/hospital_locations.rds")
hospital_locations
population <- usmap::countypop
hospital_count <- hospitals %>%
count(state, original_county)
state_hospitals <- population %>%
left_join(hospital_count, by = c("abbr" = "state", "county" = "original_county")) %>%
mutate(hospitals = ifelse(is.na(n), 0, n)) %>%
select(
fips,
state = abbr,
county,
population = pop_2015,
hospitals
)
state_hospitals
set.seed(100)
hospital_sample <- state_hospitals %>%
sample_frac(0.3)
hospital_sample
model <- lm(hospitals ~ population, data = hospital_sample)
write_rds(model, "../data/model.rds")
summary(model)
Call:
lm(formula = hospitals ~ population, data = hospital_sample)
Residuals:
Min 1Q Median 3Q Max
-10.7909 -0.7106 0.0908 0.2498 8.3805
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.027e-01 3.942e-02 17.83 <2e-16 ***
population 7.907e-06 1.182e-07 66.90 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.149 on 941 degrees of freedom
Multiple R-squared: 0.8263, Adjusted R-squared: 0.8261
F-statistic: 4475 on 1 and 941 DF, p-value: < 2.2e-16
predictions <- predict(model,
newdata = state_hospitals,
interval = "prediction") %>%
as_tibble() %>%
mutate_all(round)
predictions
hospital_results <- state_hospitals %>%
bind_cols(predictions) %>%
mutate(result = case_when(
hospitals < lwr ~ -1,
hospitals > upr ~ 1,
TRUE ~ 0)) %>%
mutate(
county_key = str_remove_all(county, " County"),
county_key = str_remove_all(county_key, " Parish"),
county_key = str_remove_all(county_key, " city"),
county_key = str_to_lower(county_key),
county_key = str_remove_all(county_key, " "),
county_key = str_replace_all(county_key, "st. ", "saint"),
)
write_rds(hospital_results, "../data/hospitals.rds")
hospital_results
State map
shapes <- read_rds("../data/shapes.rds") %>%
mutate(
county_key = str_to_lower(county),
county_key = str_remove_all(county_key, " "),
county_key = str_replace_all(county_key, "st. ", "saint")
)
shapes
county_hospitals <- hospital_results %>%
select(- county) %>%
left_join(shapes, by = c("state", "county_key")) %>%
filter(state != "PR", state != "AK") %>%
filter(!is.na(county))
write_rds(county_hospitals, "../data/county_hospitals.rds", compress = "gz")
county_hospitals
states <- county_hospitals %>%
group_by(state) %>%
summarise() %>%
pull()
under <- "#CC79A7"
over <- "#0072B2"
at_level <- "#008b00"
hospital_color <- "#F0E442"
st <- c("LA")
#st <- states
include <- c(-1, 1, 0)
add_hospitals <- 0
initial_map <- leaflet() %>%
addProviderTiles(providers$CartoDB) %>%
addLegend("bottomright",
color = c(under,over ,at_level, hospital_color),
labels = c("Less hopitals than expected",
"More hospitals than expected", "Within Range",
"Hospital Location"),
title = "Legend",opacity = 0.5)
locations <- hospital_locations %>%
filter(state %in% st) %>%
select(longitude, latitude) %>%
mutate_all(round, 3) %>%
count(longitude, latitude) %>%
mutate(popup = paste0("Hospitals: ", n))
counties <- county_hospitals %>%
filter(state %in% st,
result %in% include) %>%
mutate(popup = paste0("<b>", county, "</b>",
"<br/>Population: ", prettyNum(population, big.mark = ","),
"<br/>Hospitals: ", hospitals,
"<br/>Expected: ", fit
)) %>%
mutate(color = case_when(
result == 0 ~ at_level,
result == 1 ~ over,
result == -1 ~ under
))
county_map <- counties %>%
transpose() %>%
map(~{
county <- .x
transpose(county$data) %>%
map(~list(
data = .x$data,
color = county$color,
popup = county$popup
))
}) %>%
flatten() %>%
reduce(
~ addPolygons(.x,
lng = .y$data$long,
lat = .y$data$lat,
color = .y$color,
popup = .y$popup,
weight = 1, fillOpacity = 0.3),
.init = initial_map)
if(add_hospitals == 1) {
county_map <- county_map %>%
addCircleMarkers(
lng = locations$longitude,
lat = locations$latitude,
radius = 1.5 *locations$n,
popup = locations$popup,
fillColor = hospital_color, color = "gray",
fillOpacity = 0.7,weight = 1)
}
county_map
NA
LS0tDQp0aXRsZTogIkFjY2VzcyB0byBjYXJlIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KYGBge3J9DQpsaWJyYXJ5KHRpZHl2ZXJzZSkNCmxpYnJhcnkocmxhbmcpDQpsaWJyYXJ5KGxlYWZsZXQpDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbHMgPC0gcmVhZF9jc3YoIi4uL2RhdGEvaG9zcGl0YWxfbGlzdC5jc3YiKQ0KaG9zcGl0YWxzDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbF9sb2NhdGlvbnMgPC0gaG9zcGl0YWxzICAlPiUNCiAgc2VsZWN0KHN0YXRlLCBsb25naXR1ZGUsIGxhdGl0dWRlKQ0Kd3JpdGVfcmRzKGhvc3BpdGFsX2xvY2F0aW9ucywgIi4uL2RhdGEvaG9zcGl0YWxfbG9jYXRpb25zLnJkcyIpDQpob3NwaXRhbF9sb2NhdGlvbnMNCmBgYA0KDQpgYGB7cn0NCnBvcHVsYXRpb24gPC0gdXNtYXA6OmNvdW50eXBvcA0KDQpob3NwaXRhbF9jb3VudCA8LSBob3NwaXRhbHMgJT4lDQogIGNvdW50KHN0YXRlLCBvcmlnaW5hbF9jb3VudHkpIA0KDQoNCnN0YXRlX2hvc3BpdGFscyA8LSBwb3B1bGF0aW9uICU+JQ0KICBsZWZ0X2pvaW4oaG9zcGl0YWxfY291bnQsIGJ5ID0gYygiYWJiciIgPSAic3RhdGUiLCAiY291bnR5IiA9ICJvcmlnaW5hbF9jb3VudHkiKSkgJT4lDQogIG11dGF0ZShob3NwaXRhbHMgPSBpZmVsc2UoaXMubmEobiksIDAsIG4pKSAlPiUNCiAgc2VsZWN0KA0KICAgIGZpcHMsIA0KICAgIHN0YXRlID0gYWJiciwgDQogICAgY291bnR5LCANCiAgICBwb3B1bGF0aW9uID0gcG9wXzIwMTUsDQogICAgaG9zcGl0YWxzDQogICAgKQ0KDQpzdGF0ZV9ob3NwaXRhbHMNCmBgYA0KDQoNCg0KDQpgYGB7cn0NCnNldC5zZWVkKDEwMCkNCg0KaG9zcGl0YWxfc2FtcGxlIDwtIHN0YXRlX2hvc3BpdGFscyAlPiUNCiAgc2FtcGxlX2ZyYWMoMC4zKQ0KDQpob3NwaXRhbF9zYW1wbGUNCmBgYA0KDQpgYGB7cn0NCm1vZGVsIDwtIGxtKGhvc3BpdGFscyB+IHBvcHVsYXRpb24sIGRhdGEgPSBob3NwaXRhbF9zYW1wbGUpDQp3cml0ZV9yZHMobW9kZWwsICIuLi9kYXRhL21vZGVsLnJkcyIpDQpzdW1tYXJ5KG1vZGVsKQ0KYGBgDQoNCmBgYHtyfQ0KcHJlZGljdGlvbnMgPC0gcHJlZGljdChtb2RlbCwgDQogIG5ld2RhdGEgPSBzdGF0ZV9ob3NwaXRhbHMsIA0KICBpbnRlcnZhbCA9ICJwcmVkaWN0aW9uIikgJT4lDQogIGFzX3RpYmJsZSgpICU+JQ0KICBtdXRhdGVfYWxsKHJvdW5kKQ0KDQpwcmVkaWN0aW9ucw0KYGBgDQoNCmBgYHtyfQ0KaG9zcGl0YWxfcmVzdWx0cyA8LSBzdGF0ZV9ob3NwaXRhbHMgJT4lDQogIGJpbmRfY29scyhwcmVkaWN0aW9ucykgJT4lDQogIG11dGF0ZShyZXN1bHQgPSBjYXNlX3doZW4oDQogICAgaG9zcGl0YWxzIDwgbHdyIH4gLTEsDQogICAgaG9zcGl0YWxzID4gdXByIH4gMSwNCiAgICBUUlVFIH4gMCkpICU+JQ0KICBtdXRhdGUoDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eSwgIiBDb3VudHkiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBQYXJpc2giKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBjaXR5IiksDQogICAgY291bnR5X2tleSA9IHN0cl90b19sb3dlcihjb3VudHlfa2V5KSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiAiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlcGxhY2VfYWxsKGNvdW50eV9rZXksICJzdC4gIiwgInNhaW50IiksDQogICAgICApICANCg0Kd3JpdGVfcmRzKGhvc3BpdGFsX3Jlc3VsdHMsICIuLi9kYXRhL2hvc3BpdGFscy5yZHMiKQ0KDQpob3NwaXRhbF9yZXN1bHRzIA0KYGBgDQoNCiMjIFN0YXRlIG1hcA0KDQpgYGB7cn0NCnNoYXBlcyA8LSByZWFkX3JkcygiLi4vZGF0YS9zaGFwZXMucmRzIikgJT4lDQogIG11dGF0ZSgNCiAgICBjb3VudHlfa2V5ID0gc3RyX3RvX2xvd2VyKGNvdW50eSksDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eV9rZXksICIgIiksDQogICAgY291bnR5X2tleSA9IHN0cl9yZXBsYWNlX2FsbChjb3VudHlfa2V5LCAic3QuICIsICJzYWludCIpDQogICkgDQoNCnNoYXBlcyANCmBgYA0KYGBge3J9DQpjb3VudHlfaG9zcGl0YWxzIDwtIGhvc3BpdGFsX3Jlc3VsdHMgICU+JQ0KICBzZWxlY3QoLSBjb3VudHkpICU+JQ0KICBsZWZ0X2pvaW4oc2hhcGVzLCBieSA9IGMoInN0YXRlIiwgImNvdW50eV9rZXkiKSkgJT4lDQogIGZpbHRlcihzdGF0ZSAhPSAiUFIiLCBzdGF0ZSAhPSAiQUsiKSAlPiUNCiAgZmlsdGVyKCFpcy5uYShjb3VudHkpKQ0KDQp3cml0ZV9yZHMoY291bnR5X2hvc3BpdGFscywgIi4uL2RhdGEvY291bnR5X2hvc3BpdGFscy5yZHMiLCBjb21wcmVzcyA9ICJneiIpDQoNCmNvdW50eV9ob3NwaXRhbHMNCmBgYA0KDQpgYGB7cn0NCnN0YXRlcyA8LSBjb3VudHlfaG9zcGl0YWxzICU+JQ0KICBncm91cF9ieShzdGF0ZSkgJT4lDQogIHN1bW1hcmlzZSgpICU+JQ0KICBwdWxsKCkNCg0KYGBgDQoNCg0KYGBge3J9DQoNCnVuZGVyIDwtICIjQ0M3OUE3Ig0Kb3ZlciA8LSAiIzAwNzJCMiINCmF0X2xldmVsIDwtICIjMDA4YjAwIg0KaG9zcGl0YWxfY29sb3IgPC0gIiNGMEU0NDIiDQpzdCA8LSBjKCJMQSIpDQojc3QgPC0gc3RhdGVzDQppbmNsdWRlIDwtIGMoLTEsIDEsIDApDQphZGRfaG9zcGl0YWxzIDwtIDANCg0KaW5pdGlhbF9tYXAgPC0gbGVhZmxldCgpICU+JQ0KICBhZGRQcm92aWRlclRpbGVzKHByb3ZpZGVycyRDYXJ0b0RCKSAlPiUNCiAgYWRkTGVnZW5kKCJib3R0b21yaWdodCIsIA0KICAgICAgICAgICAgY29sb3IgID0gYyh1bmRlcixvdmVyICxhdF9sZXZlbCwgaG9zcGl0YWxfY29sb3IpLCANCiAgICAgICAgICAgIGxhYmVscyA9IGMoIkxlc3MgaG9waXRhbHMgdGhhbiBleHBlY3RlZCIsDQogICAgICAgICAgICAgICAgICAgICAgICJNb3JlIGhvc3BpdGFscyB0aGFuIGV4cGVjdGVkIiwgIldpdGhpbiBSYW5nZSIsIA0KICAgICAgICAgICAgICAgICAgICAgICAiSG9zcGl0YWwgTG9jYXRpb24iKSwgDQogICAgICAgICAgICB0aXRsZSAgPSAiTGVnZW5kIixvcGFjaXR5ID0gMC41KQ0KDQpsb2NhdGlvbnMgPC0gaG9zcGl0YWxfbG9jYXRpb25zICU+JQ0KICBmaWx0ZXIoc3RhdGUgJWluJSBzdCkgJT4lDQogIHNlbGVjdChsb25naXR1ZGUsIGxhdGl0dWRlKSAlPiUNCiAgbXV0YXRlX2FsbChyb3VuZCwgMykgJT4lDQogIGNvdW50KGxvbmdpdHVkZSwgbGF0aXR1ZGUpICU+JQ0KICBtdXRhdGUocG9wdXAgPSBwYXN0ZTAoIkhvc3BpdGFsczogIiwgbikpDQoNCmNvdW50aWVzIDwtIGNvdW50eV9ob3NwaXRhbHMgJT4lIA0KICBmaWx0ZXIoc3RhdGUgJWluJSBzdCwNCiAgICAgICAgIHJlc3VsdCAlaW4lIGluY2x1ZGUpICU+JQ0KICBtdXRhdGUocG9wdXAgPSBwYXN0ZTAoIjxiPiIsIGNvdW50eSwgIjwvYj4iLA0KICAgICAgICAgICAgICAgICAgICAgICAgIjxici8+UG9wdWxhdGlvbjogIiwgcHJldHR5TnVtKHBvcHVsYXRpb24sIGJpZy5tYXJrID0gIiwiKSwgDQogICAgICAgICAgICAgICAgICAgICAgICAiPGJyLz5Ib3NwaXRhbHM6ICIsIGhvc3BpdGFscywNCiAgICAgICAgICAgICAgICAgICAgICAgICI8YnIvPkV4cGVjdGVkOiAiLCBmaXQNCiAgICAgICAgICAgICAgICAgICAgICAgICkpICU+JSAgDQogIG11dGF0ZShjb2xvciA9IGNhc2Vfd2hlbigNCiAgICByZXN1bHQgPT0gIDAgIH4gYXRfbGV2ZWwsDQogICAgcmVzdWx0ID09ICAxICB+IG92ZXIsDQogICAgcmVzdWx0ID09IC0xICB+IHVuZGVyDQogICkpIA0KDQpjb3VudHlfbWFwIDwtIGNvdW50aWVzICU+JQ0KICB0cmFuc3Bvc2UoKSAlPiUNCiAgbWFwKH57DQogICAgY291bnR5IDwtIC54DQogICAgIHRyYW5zcG9zZShjb3VudHkkZGF0YSkgICU+JQ0KICAgICAgIG1hcCh+bGlzdCgNCiAgICAgICAgIGRhdGEgPSAueCRkYXRhLCANCiAgICAgICAgIGNvbG9yID0gY291bnR5JGNvbG9yLA0KICAgICAgICAgcG9wdXAgPSBjb3VudHkkcG9wdXANCiAgICAgICAgICkpDQogICAgfSkgJT4lDQogIGZsYXR0ZW4oKSAlPiUNCiAgcmVkdWNlKA0KICAgIH4gYWRkUG9seWdvbnMoLngsIA0KICAgICAgICAgICAgICAgICAgbG5nID0gLnkkZGF0YSRsb25nLCANCiAgICAgICAgICAgICAgICAgIGxhdCA9IC55JGRhdGEkbGF0LCANCiAgICAgICAgICAgICAgICAgIGNvbG9yID0gLnkkY29sb3IsIA0KICAgICAgICAgICAgICAgICAgcG9wdXAgPSAueSRwb3B1cCwNCiAgICAgICAgICAgICAgICAgIHdlaWdodCA9IDEsIGZpbGxPcGFjaXR5ID0gMC4zKSwgDQogICAgLmluaXQgPSBpbml0aWFsX21hcCkgIA0KDQoNCmlmKGFkZF9ob3NwaXRhbHMgPT0gMSkgew0KICBjb3VudHlfbWFwIDwtIGNvdW50eV9tYXAgJT4lDQogICAgYWRkQ2lyY2xlTWFya2VycygNCiAgICAgIGxuZyA9IGxvY2F0aW9ucyRsb25naXR1ZGUsIA0KICAgICAgbGF0ID0gbG9jYXRpb25zJGxhdGl0dWRlLA0KICAgICAgcmFkaXVzID0gMS41ICpsb2NhdGlvbnMkbiwNCiAgICAgIHBvcHVwID0gbG9jYXRpb25zJHBvcHVwLA0KICAgICAgZmlsbENvbG9yID0gaG9zcGl0YWxfY29sb3IsIGNvbG9yID0gImdyYXkiLCANCiAgICAgIGZpbGxPcGFjaXR5ID0gMC43LHdlaWdodCA9IDEpDQogIH0gDQoNCmNvdW50eV9tYXANCiAgDQpgYGANCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0K